Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add MathProgBase solver interface. #9

Merged
merged 23 commits into from
Feb 13, 2018
Merged

Conversation

tkoolen
Copy link
Collaborator

@tkoolen tkoolen commented Feb 2, 2018

I know MathOptInterface is coming soon, but I want to try this with JuMP right now.

No tests yet, but most things seem to be working. Demo code:

using OSQP, JuMP

solver = OSQPSolver()
MathProgBase.setparameters!(solver; Silent = true)
m = Model(solver = solver)
@variable(m, x)
@variable(m, y)
@constraint(m, 1 <= x <= 4)
@constraint(m, y >= 3)
@objective(m, Min, x^2 + y^2)
status = solve(m)
@show getvalue(x)
@show getvalue(y)
@show getobjectivevalue(m)

For ease of implementation, I decided to store copies of P, q, etc. in OSQPModel and call OSQP.setup! in MathProgBase.optimize!. None of the updating infrastructure is currently used. Still, it's a start.

@mlubin
Copy link
Collaborator

mlubin commented Feb 2, 2018

Looks reasonable. You should run the relevant functions in MathProgBase/test/quadprog.jl as part of the unit tests.

But I'd add that the important MOI infrastructure for QPs is now in place in JuMP, and we could use some early testers.

settings::Dict{Symbol,Any}
end

OSQPSolver() = OSQPSolver(Dict{Symbol, Any}())
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe name this OSQPMPBSolver to avoid the imminent conflict with MOI?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see this is still the common approach in most of the current MathProgBase interfaces, e.g., EcosSolverInterface.jl.

@mlubin Are you planning to change the other solvers as well or to just drop the MathProgBase interfaces altogether?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

MPB interfaces will all be moved into separate packages when we start releasing the new MOI interfaces. It's up to developer discretion whether the MOI interface will sit in the same package as the direct solver wrapper or in a separate package.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then it could be better to move this interface to a separate package called MathProgBaseOSQP.jl as done in MathProgBaseMosek.jl. In this way we could keep OSQPSolver in both MPB and MOI interfaces and avoid using both them at the same time. Is this the purpose?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems a bit silly to create a package just for 200 lines of code that are going to be deprecated soon. (I'm reconsidering if we should do this for all solvers.) If we remove the line:

import .OSQPSolverInterface: OSQPSolver

then there's nothing much to worry about with having both MPB and MOI interfaces in the same package.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Renamed the module to OSQPMathProgBaseInterface (solver is still named OSQPSolver). No longer exporting OSQPSolver, but exporting OSQPMathProgBaseInterface instead. What do you think?

@bstellato
Copy link
Collaborator

Thank you for this. I guess the next steps are to:

  • Complete the remaining functions that are compatible with OSQP
  • Implement the problem data updates properly (this gives OSQP the real speedups by the way). I think this can be done by
    • keeping track of the problem data as an internal copy of P, q, A, l, u
    • calling the osqp function updates whenever possible, e.g., when setobj!(m, c) is called we call OSQP.update!(model, q=c)
    • having a boolean flag to define whether a new setup is required when optimize! is called. This is useful when funcitons like addvar! are called.
  • Add the quadprog tests as in Gurobi.jl

@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 5, 2018

You should run the relevant functions in MathProgBase/test/quadprog.jl as part of the unit tests.

Since OSQP doesn't support quadratic constraints, neither quadprogtest nor qpdualtest is really applicable. There's also the additional issue of lack of support for variable bounds (do you agree with my handling of that by the way, or should they automatically be turned into constraints?). I opted to just copy over the relevant parts of quadprogtest for now. The dual-related functionality should still be tested somehow though.

@codecov-io
Copy link

codecov-io commented Feb 5, 2018

Codecov Report

Merging #9 into master will increase coverage by 8.22%.
The diff coverage is 89.35%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master       #9      +/-   ##
==========================================
+ Coverage   71.69%   79.91%   +8.22%     
==========================================
  Files           3        4       +1     
  Lines         272      488     +216     
==========================================
+ Hits          195      390     +195     
- Misses         77       98      +21
Impacted Files Coverage Δ
src/OSQP.jl 83.33% <ø> (ø) ⬆️
src/interface.jl 69.54% <ø> (+0.82%) ⬆️
src/mpbinterface.jl 89.35% <89.35%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update eb1ea2d...5094552. Read the comment docs.

sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf,solver)
@test sol.status == :Optimal
@test isapprox(sol.objval, 130/70, atol=1e-6)
@test isapprox(norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), 0.0, atol=1e-6)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm surprised the test passes with these tolerances.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oops, missed the setting for eps_abs above.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, polish = true is more essential. Just setting eps_abs to even 1e-8 results in test failure.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Polish shouldn't be necessary, but you may need to set eps_rel (osqp/osqp#40) to a tiny value.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you take the relative tolerance into account it might work without polishing:

@test isapprox(norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), rtol=1e-6, atol=1e-6)

REQUIRE Outdated
@@ -1,3 +1,4 @@
BinDeps
julia 0.6
Compat
MathProgBase 0.7
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's always a good idea to put an upper bound on the MPB requirement.

@mlubin
Copy link
Collaborator

mlubin commented Feb 5, 2018

There's also the additional issue of lack of support for variable bounds (do you agree with my handling of that by the way, or should they automatically be turned into constraints?).

Depends on whether you want to call OSQP on JuMP models that have variables with bounds.

@bstellato
Copy link
Collaborator

Since OSQP doesn't support quadratic constraints, neither quadprogtest nor qpdualtest is really applicable. There's also the additional issue of lack of support for variable bounds (do you agree with my handling of that by the way, or should they automatically be turned into constraints?). I opted to just copy over the relevant parts of quadprogtest for now. The dual-related functionality should still be tested somehow though.

I think people use variable bounds a lot. We could keep track of them in the problem data as

    lb <= Ax <= ub
    l <= x <= u

Then, when we call optimize! we stack them together. I can add it later today or tomorrow.

Only allow :Min as the sense for now, as the current implementation can lead to
confusing/erroneous behavior.

Add freemodel!, copy, getsense, setsense!.
@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 5, 2018

Re: #9 (comment): I started doing that initially, but I got a little worried that the interpretation of the dual variables would become confusing/nonstandard.

We also need better handling of sense. For now, I've restricted sense to be :Min in my latest commit in order to at least avoid blatantly wrong behavior. I'm not sure what the JuMP/MathProgBase convention is for handling duals when sense == :Max, so guidance would be appreciated. Also, the fact that loadproblem! contains the sense but not the 'P' matrix would have made subsequent calls to setquadobj! do the wrong thing for sense == :Max.

@mlubin
Copy link
Collaborator

mlubin commented Feb 5, 2018

For now, I've restricted sense to be :Min in my latest commit in order to at least avoid blatantly wrong behavior. I'm not sure what the JuMP/MathProgBase convention is for handling duals when sense == :Max, so guidance would be appreciated.

JuMP has lots of tests for duals on LPs with both min/max sense. Not sure if there are specific tests for QPs. You should really feel free to only implement the parts that you need, because all of this is going to be tossed out after JuMP 0.18. Rejecting Max problems is reasonable.

@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 6, 2018

Improved test coverage. Also discovered and fixed a minor problem with OSQP.warm_start! along the way.

@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 6, 2018

Don't merge yet; I suspect there's something wrong with the quadratic objective methods.

@bstellato
Copy link
Collaborator

Thanks for adding these tests. Tomorrow I will merge my changes as well to complete the interface.

@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 7, 2018

Fixed the objective bug and added yet more tests that displayed the problem. I'm successfully using OSQP through JuMP for my intended application now.

end
end

function Base.resize!(model::OSQPModel, n, m)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the resize! function really necessary? it seems that it does not properly resize the problem but it just creates zero matrices P and A and resizes the vectors. I guess this is needed in the loadproblem! function where you call copy!. Is this necessary/more efficient than just using copy or deepcopy for the data?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FTR: answered on gitter. Yes, efficiency.

@tkoolen tkoolen mentioned this pull request Feb 9, 2018
22 tasks
Copy link
Collaborator Author

@tkoolen tkoolen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In general, there will be quite a lot of allocations when using the additional functionality. You could add some fields to the model to avoid those allocations. Or (and that's absolutely fine) don't bother and we'll handle it properly in the MathOptInterface version.

test/runtests.jl Outdated
# "primal_infeasibility.jl",
# "unconstrained.jl",
# "warm_start.jl",
# "update_matrices.jl",
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reenable

module OSQPMathProgBaseInterface

using OSQP: Model, Results, setup!, solve!, update!, clean!, update_settings!, warm_start!
importall MathProgBase.SolverInterface
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

end
end

function get_qp_variables(model::OSQPMathProgModel)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not just have this be numvar?

return n
end

get_qp_constraints(model::OSQPMathProgModel) = size(model.A, 1)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

numconstr?

checksolved(model::OSQPMathProgModel) = isdefined(model, :results) || error("Model has not been solved.")

# Reset problem when setup has to be performed again
resetproblem(model::OSQPMathProgModel) = (model.results = nothing; model.perform_setup=true)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This shouldn't work. model.results can't store nothing.


# Negate cost if necessary
if model.sense == :Max
model.q *= -1
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Allocation could be avoided, but not a huge deal.

getconstrmatrix(model::OSQPMathProgModel) = model.A


function addvar!(model::OSQPMathProgModel, constridx, constrcoef, l, u, objcoef)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Lots of allocations, to the point that it may not be beneficial to support this. I don't think JuMP even uses this.

@boundscheck length(colidx) == nterms || error()

# Check if only the values have changed
if isdefined(model, :P)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is always true.

# Check if only the values have changed
if isdefined(model, :P)
Pi, Pj, Px = findnz(model.P)
if (rowidx == Pi) & (colidx == Pj) & !model.perform_setup
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that ordering of rowidx and colidx shouldn't matter and

Duplicate index sets (i,j) are accepted and will be summed together

this is kind of a crude check.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have left it for now since I cannot find an easier way to check if the sparsity pattern of P stays the same. If you have any suggestions feel free to update it.


# Change sign if maximizing
if model.sense == :Max
model.P *= -1
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Expensive way of doing this.

@bstellato
Copy link
Collaborator

Damn, I pushed before finishing some stuff!!! Sorry. I had to rush out to finish some other things. I wanted to rebase everything before pushing it. I will adjust those things tomorrow ok? (Have lots of IKEA stuff to put together...).

Thanks anyway for the performance tips. I am not totally familiar with Julia performance guidelines yet. I haven't used it that much until very recently.

@bstellato
Copy link
Collaborator

I have fixed the interface. It should now support incremental model building, infeasibility certificates, Max or Min problems and problem data updates. There was a problem with the dual variables sign before. I have removed a couple of unnecessary allocations but I am quite sure we could do something better memory-wise. However, we might just leave it as it is and focus on performance improvements for the MathOptInterface. Let me know what you think.

@tkoolen
Copy link
Collaborator Author

tkoolen commented Feb 13, 2018

Looks good to me!

@bstellato bstellato merged commit 69d25e7 into osqp:master Feb 13, 2018
@tkoolen tkoolen deleted the tk/mathprogbase branch February 13, 2018 03:16
bstellato pushed a commit that referenced this pull request May 17, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants